Modelos lineales

Regresión y factoriales (ANOVA y ANCOVA)

Introducción

  • Definición

  • Tipos

  • Supuestos

  • Flujo de trabajo

Definición

Relación lineal entre dos o más variables, dependiendo una (variable dependiente) de la(s) otra(s) (variable(s) explicativa(s))

Modelo matemático lineal: y = 2x o y = 3 - 2x

Objetivos

  1. Descripción de la relación

  2. Inferencia, o generalización de los resultados

  3. Predicciones

Tipos

Variable dependiente: siempre cuantitativa

Variables explicativas:

  • Cuantitativas: modelos de regresión (simple o múltiple)

  • Cualitativas o categóricas: ANOVA

  • Cuantitativas y categóricas: ANCOVA

Supuestos del modelo

  1. Independencia

  2. Linealidad

  3. Normalidad

  4. Homocedasticidad

Flujo de trabajo

Pasos a seguir análisis de datos con modelos lineales (Cayuela y de la Cruz, 2022)

Formular problema correctamente

  1. Comprender problema y contexto

  2. Comprender objetivo del estudio

  3. Plantear problema en términos estadísticos

  4. Entender bien los datos

Regresión simple y regresión múltiple

Estructura

Ecuación regresión

β0: ordenada en el origen

β1: pendiente

Ɛi: residuos

Ajuste modelo lineal

Función lm(y ~ x)

Ejemplo

Pregunta: Masa corporal especies pingüinos y cómo esta es influenciada por la longitud de la aleta de los pingüinos

Hipótesis

H0 (hipótesis nula): no existe relación entre la masa corporal de los pingüinos y la longitud de sus aletas

H0 : β1 = 0

H1 (hipótesis alternativa): pendienta de la recta que relaciona masa corporal con longitud aleta es distinta de cero

H0 : β1 ≠ 0

1. Leer los datos en R y explorar gráficamente la relación entre las variables

library(palmerpenguins)
data(penguins)
head(penguins)
# A tibble: 6 × 8
  species island    bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
  <fct>   <fct>              <dbl>         <dbl>             <int>       <int>
1 Adelie  Torgersen           39.1          18.7               181        3750
2 Adelie  Torgersen           39.5          17.4               186        3800
3 Adelie  Torgersen           40.3          18                 195        3250
4 Adelie  Torgersen           NA            NA                  NA          NA
5 Adelie  Torgersen           36.7          19.3               193        3450
6 Adelie  Torgersen           39.3          20.6               190        3650
# ℹ 2 more variables: sex <fct>, year <int>
library(ggplot2)
ggplot(data = penguins, aes(x = flipper_length_mm, y = body_mass_g)) +  
  geom_point(color = "blue", size = 3) +      
  labs(x = "Longitud aleta (mm)",         
       y = "Masa corporal (g)",                   
      title = "Relación entre masa corporal y longitud aleta "   
  ) +
  theme_minimal()

2. Ajuste modelo de regresión

lm_reg1 <- lm(body_mass_g ~ flipper_length_mm, data = penguins)

3. Variación explicada por el modelo y resolución hipótesis

Tabla del análisis de la varianza
anova(lm_reg1)
Analysis of Variance Table

Response: body_mass_g
                   Df    Sum Sq   Mean Sq F value    Pr(>F)    
flipper_length_mm   1 166452902 166452902  1070.7 < 2.2e-16 ***
Residuals         340  52854796    155455                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Suma de cuadrados: variabilidad de la variable y

Coeficiente de determinación (R2): bondad de ajuste del modelo (0-1)

4. Interpretación del modelo

summary(lm_reg1)

Call:
lm(formula = body_mass_g ~ flipper_length_mm, data = penguins)

Residuals:
     Min       1Q   Median       3Q      Max 
-1058.80  -259.27   -26.88   247.33  1288.69 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)       -5780.831    305.815  -18.90   <2e-16 ***
flipper_length_mm    49.686      1.518   32.72   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 394.3 on 340 degrees of freedom
  (2 observations deleted due to missingness)
Multiple R-squared:  0.759, Adjusted R-squared:  0.7583 
F-statistic:  1071 on 1 and 340 DF,  p-value: < 2.2e-16

y = 49.686x - 5780.831

Fiabilidad estimadores para predecir

R2 y R2 ajustado

5. Predicciones con el modelo

library(broom)
library(tidyverse)
library(dplyr)
preds <- augment(lm_reg1, type.predict = "response", se_fit = TRUE)
preds <- preds %>%
  mutate(lower = .fitted - 1.96 * .se.fit,upper = .fitted + 1.96 * .se.fit)
  
ggplot(data = preds, aes(x = flipper_length_mm, y = body_mass_g)) +  
  geom_point(color = "blue", size = 2) +     
  geom_line(aes(y = .fitted)) +
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2) +
  labs(x = "Longitud aleta (mm)",         
       y = "Masa corporal (g)",                   
      title = "Predicciones modelo de regresión"   
  ) +
  theme_minimal()

6. Revisión supuestos del modelo

Supuestos de linealidad, normalidad y homocedasticidad

par(mfrow=c(2,2))
plot(lm_reg1)

Normalidad:

shapiro.test(residuals(lm_reg1))

    Shapiro-Wilk normality test

data:  residuals(lm_reg1)
W = 0.99301, p-value = 0.1123

Linealidad:

library(lmtest)
resettest(lm_reg1)

    RESET test

data:  lm_reg1
RESET = 15.998, df1 = 2, df2 = 338, p-value = 2.3e-07

Homocedasticidad:

library(car)
ncvTest(lm_reg1)
Non-constant Variance Score Test 
Variance formula: ~ fitted.values 
Chisquare = 2.177953, Df = 1, p = 0.14

No se cumplen supuestos del modelo: transformación variables (logarítmica, raíz cuadrada, Box-Cox, log-log, etc.)

Regresión múltiple

yi = β0 + β1x1 + … + β1xn + Ɛi

Ɛi ~N(0,σ2)

Suma de cuadrados tipo I o III

Suma de cuadrados: tipos (Cayuela y de la Cruz, 2022

Modelos factoriales: ANOVA y ANCOVA

Análisis de la varianza (ANOVA)

Hipótesis: comparación medias distintos grupso (niveles de factor)

H0: medias de los grupos son iguales

H0: µa = µb = … = µ

H1: al menos una media es distinta a las demás

H0: Ǝµj ≠ µ

Ejemplo

Pregunta: Diferencias longitud pico tres especies de pingüino

Hipótesis

H0 (hipótesis nula): no existe diferencias en las medias de la longitud de los picos para las tres especies

H0 : µadelie = µgentoo = µchinstrap = µ

H1 (hipótesis alternativa): la media de al menos una de las tres especies es distinta a las demás

H0 : Ǝµj ≠ µ

1. Explorar gráficamente la relación entre las variables

ggplot(data = penguins, aes(x = species, y = bill_length_mm, fill = species)) +  
  geom_boxplot(alpha = 0.8, color = "gray30") +
  scale_fill_manual(values = c("#F8BBD0", "#FFF9C4", "#B3E5FC")) +
  labs(x = "Species",         
       y = "Longitud pico (mm)",                   
      title = "Diferencias longitud pico entre especies"   
  ) +
  theme_minimal()

2. Ajuste modelo factorial

penguins$species <- as.factor(penguins$species)
lm_reg2 <- lm(bill_length_mm ~ species, data = penguins)

3. Variación explicada por el modelo y resolución hipótesis

Tabla de análisis de la varianza
anova(lm_reg2)
Analysis of Variance Table

Response: bill_length_mm
           Df Sum Sq Mean Sq F value    Pr(>F)    
species     2 7194.3  3597.2   410.6 < 2.2e-16 ***
Residuals 339 2969.9     8.8                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

4. Interpretación del modelo y comparaciones post-hoc

summary(lm_reg2)

Call:
lm(formula = bill_length_mm ~ species, data = penguins)

Residuals:
    Min      1Q  Median      3Q     Max 
-7.9338 -2.2049  0.0086  2.0662 12.0951 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)       38.7914     0.2409  161.05   <2e-16 ***
speciesChinstrap  10.0424     0.4323   23.23   <2e-16 ***
speciesGentoo      8.7135     0.3595   24.24   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.96 on 339 degrees of freedom
  (2 observations deleted due to missingness)
Multiple R-squared:  0.7078,    Adjusted R-squared:  0.7061 
F-statistic: 410.6 on 2 and 339 DF,  p-value: < 2.2e-16

Variables indicadoras o variables dummy

yi,j = β1x1 + β2x2 + β3x3 + Ɛi

Un nivel es la coordenada en el origen (nivel de referencia)

 contrasts(penguins$species)
          Chinstrap Gentoo
Adelie            0      0
Chinstrap         1      0
Gentoo            0      1

Ecuación: yi = 38.79 + 10.04xi + 8.71x2

library(multcomp)
posthoc_species <- glht(lm_reg2, linfct=mcp(species="Tukey"))
summary(posthoc_species)

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: lm(formula = bill_length_mm ~ species, data = penguins)

Linear Hypotheses:
                        Estimate Std. Error t value Pr(>|t|)    
Chinstrap - Adelie == 0  10.0424     0.4323  23.232  < 0.001 ***
Gentoo - Adelie == 0      8.7135     0.3595  24.237  < 0.001 ***
Gentoo - Chinstrap == 0  -1.3289     0.4473  -2.971  0.00872 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)

5. Representación gráfica diferencias significativas

ggplot(data = penguins, aes(x = species, y = bill_length_mm, fill = species, color = species)) +  
  geom_boxplot(alpha = 0.8, color = "black") +
  geom_jitter(width = 0.15, alpha = 0.6, size = 2, shape=21, color="black") +
  scale_fill_manual(values = c("#F8BBD0", "#FFF9C4", "#B3E5FC")) +
  scale_color_manual(values = c("#F8BBD0", "#FFF9C4", "#B3E5FC")) +
 labs(x = "Species",         
       y = "Longitud pico (mm)",                   
      title = "Diferencias longitud pico entre especies"   
  ) +
  theme_minimal() +
  geom_text(data = penguins[2,], aes(y = 65, label = "a"), color="black", 
            size=5, nudge_x = 0) +
  geom_text(data = penguins[2,], aes(y = 65, label = "b"), color="black", 
            size=5, nudge_x = 1) +
  geom_text(data = penguins[2,], aes(y = 65, label = "c"), color="black", 
            size=5, nudge_x = 2)

6. Revisión supuestos del modelo

par(mfrow=c(2,2))
plot(lm_reg2)

Análisis de la covarianza (ANCOVA)

Hipótesis: homogeneidad de pendientes

H0: las pendientes son iguales

H1: al menos una pendiente es distinta a las demás